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Abstract. This article investigates entanglement of the motional states of massive 
coupled oscillators. The specific realization of an idealized diatomic molecule in one- 
dimension is considered, but the techniques developed apply to any massive particles 
with two degrees of freedom and a quadratic Hamiltonian. We present two methods, 
one analytic and one approximate, to calculate the interatomic entanglement for 
Gaussian and non-Gaussian pure states as measured by the purity of the reduced 
density matrix. The cases of free and trapped molecules and hetero- and homonuclear 
molecules are treated. In general, when the trap frequency and the molecular frequency 
are very different, and when the atomic masses are equal, the atoms are highly- 
entangled for molecular coherent states and number states. Surprisingly, while the 
interatomic entanglement can be quite large even for molecular coherent states, the 
covariance of atomic position and momentum observables can be entirely explained by 
a classical model with appropriately chosen statistical uncertainty. 
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1. Introduction 

Harmonically-coupled massive oscillators provide a suitable model for many physical 
systems that are employed or proposed for quantum information processing with 
continuous variables. For example, both the longitudinal modes [1] and the transverse 
modes [2] of ions in Paul traps can be treated (in certain regimes, at least ap- 
proximately) as coupled oscillators. Interactions between oscillators, either direct or 
mediated by external elements or fields, lead to entanglement between the oscillators 
that in principle could be externally controlled or extracted. In particular, the 
characterization of entanglement in harmonic chains [3], one-dimensional arrays of 
coupled harmonic oscillators, is a paradigm that has attracted sustained attention. 
Beside quantum information processing, an additional motivation for investigating 
systems with harmonic lattice Hamiltonians is the study of the role of entanglement 
in phase transitions 

Most research on massive coupled oscillators has employed Gaussian states, the 
workhorse of continuous variable quantum information theory. This extensive use of 
Gaussian states is physically motivated: many dynamical processes result in Gaussian 
states. In particular, the ground states and thermal states of massive coupled oscillators 
with quadratic Hamiltonians are exactly Gaussian. Their use is also theoretically 
convenient because many results for Gaussian states in quantum optics can be carried 
over to massive oscillators. More generally, there exists a correspondence (up to 
local coordinate transformations) between covariance matrices and Gaussian states 
that allows for deep mathematical analysis, including complete characterization of 
bipartite entanglement and partial characterization of multipartite entanglement O [6] . 
In contrast, one goal of this paper to to provide methods for studying the entanglement 
of non- Gaussian continuous variable states. Advances in theory and experiment have 
led to an increased interest in quantum information processing with non-Gaussian states, 
and this work investigates the simplest special case of two-mode pure states. 

Our results apply directly to any case of massive particles with two degrees of 
freedom and a Hamiltonian quadratic in position, but for conceptual ease and clarity 
we structure our discussion in the language of a diatomic molecule: two distinguishable 
atoms interacting via a quadratic potential in one-dimension. This can be thought of 
as an example of the shortest harmonic chain, but the analysis will not be restricted 
to the case of equal masses as is usually considered in such systems. Two different 
Hamiltonians are considered simultaneously: a molecule trapped in a harmonic potential 
and an untrapped molecule with a Gaussian wave packet. The properties of two- 
mode Gaussian states are well-known from a variety of contexts in quantum optics 
and continuous- variable quantum information theory [3 [H [9], HUl [H] . For the ground 
state and coherent states, we will translate our model into the standard language of 
covariance matrices and logarithmic negativity. However, for non-Gaussian pure states, 
the entanglement will be quantified in terms of the purity of the reduced density matrix. 
Two methods, one more suitable for analytic calculations and one more suitable for 
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numerical simulations, will be provided. 

Entanglement is a notion that depends on the observables one uses to describe a 
system [12], and the coupled oscillator model allows this connection to be explored in 
an exactly solvable model. Because of the coupling interaction, the "normal" molecular 
observables are more convenient for studying and controlling the system properties and 
dynamics than the the "natural" atomic observables. However, one could imagine that 
atomic observables like position and momentum are still physically accessible, perhaps 
through some independent coupling with internal atomic structure. Assuming that both 
atomic observables and molecular observables can form a complete set of operationally- 
accessible interactions and measurements, one can talk about entanglement with respect 
to the tensor product structure induced by either the molecular set or the atomic 
set of observables (see also the discussion in Ref. pH]). Additionally, by studying 
the local unitary operators acting on the atomic tensor product structure, one can 
identify equivalence classes of equally-entangled states and classes of Hamiltonians that 
lead to entanglement-equivalent dynamics. For example, these methods demonstrate 
immediately that coherent states and the ground state have the same entanglement 
between atoms, and that the dynamics of the entanglement is independent of the linear 
terms in the Hamiltonian. 

We can also interpret our results for the entanglement of atoms in a diatomic 
molecule as an example of continuous- variable entanglement constrained by conservation 
laws. In this perspective, we find that the entanglement between massive oscillators 
has two determining factors, one dynamic and one kinematic in origin. First, the 
Hamiltonian is diagonal in the molecular observables, so the center-of-mass and relative 
modes are not mixed by the dynamics. As a result, the ground state, coherent states, 
and number states are separable in the molecular observables. Some dynamic parameter, 
e.g. the ratio of the molecular frequency to the trap frequency, will therefore set one scale 
in the analytic formulas for entanglement. Additionally, the transformation between 
atomic and molecular coordinates induces a kind of purely kinematic squeezing of 
the two-particle wave function. The motional entanglement created by this kind of 
wave packet squeezing is mathematically similar to phase space squeezing in quantum 
optics. This effect was perhaps first noted by Fan and Klauder [Zj who, inspired 
by the original EPR paper, studied eigenstates of the relative motion of two-particle 
systems and constructed two mode entangled states in analogy to photonic two mode 
squeezed states [14j. Subsequent work by Fan generahzed these considerations to the 
entangled state representation for two unequal masses [151 US] ■ A similar entanglement 
mechanism occurs in the reflected modes of two-particle scattering systems [171 HH] and 
the same mechanism can be found for wave packet entanglement in photoionization [19] , 
spontaneous emission [20], and other disassociation processes [21]. 

In terms of diatomic systems, possible physical realizations could include two ions 
in a linear Paul trap or cold polar diatomic molecules in an optical trap. For example, 
entanglement swapping between internal atomic degrees of freedom and molecular 
degrees of freedom has already been demonstrated for two pairs of oscillating ions in a 
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linear Paul trap [22] and novel schemes for entangling transverse modes ion traps have 
been proposed [23]. A straightforward implementation scheme for measuring this kind 
of entanglement would require independent access and measurement of both molecular 
and atomic canonical observables (i.e., position and/or momentum). Alternatively, 
if the molecule could be disassociated with a strong pulse that does not change the 
original spatial distribution of the wave packet very much, the ratio of the width of the 
conditional wave packet to the single particle wave packet takes exactly the same value 
as the purity of the reduced density matrix [191 EI]- However, either approach would 
require measurement resolution finer than the scale of the wave packet variation. 

Unfortunately, as discussed in the conclusion, even if two-particle spatial covariance 
measurements were accurate enough to quantify the interatomic entanglement, such 
measurements cannot establish the 'quantumness' of the correlations. The correlations 
between atoms revealed by spatial measurements does not exceed those that are 
possible in some classical system with statistical correlations. To overcome this, 
several schemes for developing Bell-type inequalities have been proposed for detecting 
entanglement in continuous variable systems, such as displaced parity operators [21] 
and pseudospin operators [25] . While these are useful theoretical discriminators of non- 
classical correlations, these schemes would appear to require full state tomography to 
reconstruct the correlations between arbitrary two- mode states [25]. Disassociation-time 
entanglement has also been proposed to measure motional entanglement of two atoms 
disassociated from a diatomic molecule p6] , but there the measured entanglement would 
be created by the disassociation pulse, and not the initial entanglement of the bound 
state, which is what interests us here. 

An additional physical motivation for this work is to study entanglement in bound 
states of strongly-interacting particles. For example, because the Moshinsky atom [27], 
a coupled oscillator model for two-electron atoms, is analytically solvable, it has been 
used to test and explore approximation schemes for multi-electron settings like Hartree- 
Fock [271 I2H1 [29] and density functional theory [SOj [311 [32]. The amount of spatial 
entanglement between electrons in the Moshinsky atom (and also in the related Hooke's 
atom model) has been shown to be a good proxy for the deviations in energy prediction 
entailed by the separability assumptions used in Hartee-Fock [33l |34l |35] and density 
functional theory [361 [ST]. The results presented in this article reproduce the direct 
calculations of the spatial entanglement of the ground state found in Ref. [37] and of 
the lowest energy eigenstates in Ref. [35] (when adjusted to account for the difficulties 
associated with entanglement of identical fermions). Our results also explain an 
entanglement symmetry first noticed in Ref. [35]: when the center-of-mass and relative 
energy scales are interchanged, the entanglement is invariant. 

The structure of this article is as follows. First, the model for the diatomic molecule 
is introduced, and some relevant properties of molecular and atomic entanglement are 
established. Then the atomic entanglement of coherent molecular states, including 
the ground state, is calculated and compared to known results for Gaussian states 
established using covariance matrix methods. In the next section, one method for 
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generating the exact expression for the atomic entanglement in an arbitrary energy 
eigenstatc is presented, and an alternate, approximate procedure is described in 
Appendix A. The final section discusses the interpretation of the correlations implied 
by atomic entanglement and how the covariance in this quantum system compares to a 
classical, statistical description of the same system. 



2. Idealized Diatomic Molecule 



The Hamiltonian for the trapped diatomic molecule can be expressed in terms of 
the "normal" or molecular center-of-mass observables {X, P} and relative observables 

{R,Q] 

where M = mi + 1712 is the total mass, = mim2/M is the reduced mass. 

These observables are related to the "natural" atomic observables by the symplectic 
transformation 

X = /xiXi + /X2X2, R = Xi-X2-i, 

P =Pl + p2, Q = fl2Pl - f^lP2, (2) 

where fii = rrii/M are the mass fractions and i is the equilibrium length of the molecule. 
The equilibrium length (or the coefficient of any linear term in the Hamiltonian) can be 
set equal to zero without changing entanglement because such a redefinition corresponds 
to a local unitary transformation in either molecular or atomic coordinates, as discussed 
below. Written in the atomic observables with i — 0, the Hamiltonian becomes 

H = -^P2 + 7^P2 + ViXi, X2), where (3) 

V = ^MQ^ + iiiii29^X^ + + iiiii29^X^ + AtiAt2(l - g'^)XiX2) ■ 

In the last fine we have expressed the ratio of the molecular frequency to the trap 
frequency as g — uu/D, and we will interpret all results for the bound molecule in terms 
of this dynamical scale parameter. 

Looking at the Hamiltonian in the two different coordinate systems, the benefit 
of the molecular observables is clear (and well-known). Cast into the language of 
entanglement, one can say the Hamiltonian is a separable operator in the molecular 
Hilbert space partition = H^^'Hc induced by the relative (denoted "r" ) and center- 
of-mass observables ( "c" ) 



H = H,.(g)Ic + Ir® He 

= hiu{ad^ + 1/2) Ic + Ip (8 hn{bb^ + 1/2). (4) 

In the last line the Hamiltonian has been written in terms of the ladder operators for 
the molecule oscillations {d,a^} and for the trap oscillations {6, 6^}: 

7 a ^ a 
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where T = a/ MVt/h and 7 = ^/Jujjjh = ^JJIifpigV are proportional to the momentum 
uncertainties of the ground state. In contrast, the Hamiltonian expressed in atomic 
observables ([3]) does not separate with respect to the atomic tensor product structure 
% = T-ii ® unless 5^ = 1 (or the unphysical case of infinite mass imbalance, /zi or 
/X2 — )■ 0). When = 1, the molecular and trap frequencies are the same and the term 
proportional to X1X2 vanishes (although 7 and F, which depend on the masses, may 
still be different). 

The energy eigenstates |m, n) = |m) ® \n) have the standard harmonic oscillator 
wave functions $m^„(r, x) = 0m('")0n(a;) = {r,x\m,n) when expanded on the spectrum 
of {R,X}: 

^ (2-m!)-i/2/f„(7r)e-^^'^^/2 



[X) 



TX 

^2\ 1/4 



(2"n!)-i/2^„(ra;)e-^ " /2 (g) 



il 

TT 

The energy eigenstates are separable, and therefore unentangled, with respect to the 
molecular tensor product structure Hr^'Hc, although out of these basis vectors one can 
construct combinations that are entangled in the molecular tensor product structure, 
e.g. \m,n) ± \n,m). One can also define two-mode molecular coherent states 

\a,(3) = D{a,a)D{b,(3)\0,0) 

^ g-|.|V2-|/3|V2 4i=|m,n) (7) 

m,n=0 * 

which have well-known physical interpretations as the "most classical" harmonic 
oscillator states (see, for example |38j) and will be useful for subsequent calculations. 
The complex number a is the displacement of the relative d mode and /3 is the 
displacement of the center-of-mass b mode. 

Any state |$) can also be represented by wave functions on the spectrum of the 
atomic position observables: 

$(xi, X2) = {Xi, X2|$) = $(xi - X2, lllXi + II2X2) (8) 

where we have used the "curly ket" notation |a;i,a;2} to indicate that these are 
generalized eigenvectors of the atomic position observables, as opposed to the molecular 
position observable eigenkets \r,x). For almost all values of kinematic and dynamic 
parameters, and for almost all quantum numbers and superpositions, the wave function 
is not separable, i.e. $(a;i,X2) 7^ 4'i{xi)(t)2{x2) for any (pi. This entanglement is evident 
from the contour plots of probabilities densities in xi,X2-space depicted in Figures 1 
and 2. Only the top two contour plots of Figure 1 depict separable states, which can be 
recognized because all marginal probabilities for a given value of one coordinate (say, 
Xi) take the same functional form in the other coordinate. Or more qualitatively, the 
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Figure 1. These contour plots depict the atomic position probability densities 
|$oo(a;i, a;2)P for four combinations of values for g and fii. Positions are measured 
in units of F. Note that /U2 = 1 — yUi- 



"principle axes" of the probability densities line up with the coordinate axes when the 
function is separable. 

To calculate the entanglement between atoms for arbitrary pure states, we will use 
the purity of the reduced matrix element P (or purity, for short): 

P = Tt^pI = J dxi{xi\pl\xi} 

— J dxidx[{xi\pi\x[}{x[\pi\xi} 

= J dxidx[dx2dx'2^{xi,X2)^*{x[,X2)^{x[,x'2)^*{xi,x'2) (9) 

where 



pi = Tr2(|$)($|) 



Entanglement in 



Massive Coupled Oscillators 



8 



4h 



fii = 1/2 



fii = 1/10 



2 - 



-4L 



4h 






Xi 





Pm = 0.82 



01 

-4 -2 



II 



2 - 



-2 - 



-4 b 




Poi = 0.384 

4-2 2 



2 - 



-2 - 



-4 L 




Poi =0.512 

-2 2 4 

Xi 



Figure 2. These contour plots depict the atomic position probability densities 
|<l>io(a;i, a;2)P for four combinations of values for g and fii. Positions are measured 
in units of F. Note that /U2 = 1 — yUi- 



dxidx'idx2^{xi, X2)^*{xi, X2)\xi}{x'i\. 



(10) 



The purity of the reduced density matrix is an entanglement measure for pure states 
and hes in the range (0, 1], with lower values meaning more entanglement (some prefer 
to use the linear entropy L = 1 — P for this reason). For comparison, in d-level discrete 
systems the purity is bounded from below by d^^. 

As with all reasonable entanglement measures on pure states, the value of P should 
not change when the state is transformed by a unitary transformation that is separable 
with respect to tensor product structure. In particular, operators that are exponentials 
of linear combinations of atomic or molecular position and momentum observables, such 
as 



(11) 
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are separable with respect to both the atomic tensor product structure Tii ® 'H2 and 
molecular tensor product structure Tir (8) 'He- Therefore translations of either the 
molecular or atomic coordinate systems do not affect the amount of entanglement. This 
explains why the equilibrium length of the molecule can be set to £ = with out changing 
any conclusions about entanglement. The time evolution operator on energy eigenstates 
is just a phase, so one can also see that stationary states have constant entanglement 
in time according to any tensor product structure. Finally, Fourier transforms are also 
local unitary operators. Therefore, one can work in atomic momentum space or atomic 
position space and calculate the same value for the entanglement between the atoms. 
We choose to work in position space. 

In the limit — 0, the term in the Hamiltonian ([1]) that leads to the center-of-mass 
trapping vanishes and our system becomes an untrapped diatomic molecule. We can 
still consider wave functions (E]) of the form ^rno{r,x) = (f)mir)4>o{x), but these are no 
longer energy eigenstates. The constant F now plays the role of an initial condition, 
not a dynamical parameter as in the trapped case. The quantity hT/\/2 can now be 
interpreted as the momentum uncertainty of the center-of-mass Gaussian wave packet at 
the moment in time (say t = 0) when the wave packet satisfies the minimum uncertainty 
relation AxAp = h/2. Of course, the center-of-mass wave packet will spread as a 
function of time 



where we use a unitless rescaled time r = T'^ht/M. Translations in center-of-mass 
position and momentum (or equivalently, changes in reference frame) do not affect the 
entanglement for the unbound molecule, as can be seen from the preceding argument. 
Therefore, the entanglement of the mth molecular vibrational state with a Gaussian 
center-of-mass momentum distribution can be calculated from $(t)m„(xi, X2) = (pmi^i — 
X2)(l>u{lJ'iXi + fi2X2,t). Note that since these states are no longer energy eigenstates, 
atomic entanglement will not be constant. In the next section we show that the 
entanglement for the unbound molecule increases monotonically as the wave packet 
spreads. 

3. Entanglement of Coherent States and the Ground State 

The ground state and coherent states of the diatomic Hamiltonian are Gaussian states 
in either the molecular or the atomic basis, and as such their entanglement properties 
can be specified by the corresponding covariance matrix. The symplectic eigenvalues of 
the partial transpose of a covariance matrix provide the separability criterion and can be 
used to calculate several measures of entanglement, such as logarithmic negativity [HIE]. 
At the end of this section, we will apply these methods. However, since our eventual 
goal is to evaluate ([9]) for general, non-Gaussian states, it will be instructive to first 
consider the purity Pajs of molecular (two-mode) coherent states \a,/3). 




(12) 
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In terms of the position coordinates $0/3 (r,x) 
Gaussian wave functions 



{r,x\a,f3), these states have 



<l>a^{r,x) = ( — 1 



1/2 



where 



ar 



^7 
1 

7^ 



TT / 



,a + a 



e ■^('^'■ag+l3^(5p)^-^{r-arr-i^(x-l3^y^jr{rag+xl3p) j^j^g-j 

ihj 



an 



V2 



{a — a*), 



+ and/3p 



ihr 

'7f 



(/3 -/?*). 



(14) 



Transforming to the particle coordinates via the symplectic transformation (|2]), the new 
wave function ^ai3{xi,X2) = ^api^i — X2, fJ^iXi + 112X2) is still Gaussian. The integral 
([9]) can be rewritten in the form 



0/3 



d ze 



-zTAz+B^z+C 



(15) 



with z = (xi, x'l, 0:2, Xg) ^ , d 7. = dxidx[dx2dx'2, y = l/2(— 7 + F /iiyU2) and 



A 



y 

y 





/rV? + 7' 

rV? + 7' 
y y 
\ y y 

I 7V + rVi/3x \ 
72a, + rVi/3x 
-72a, + rV2/3x 
\ -72a, + rV2/3x / 

2(7V^ + r^/3.^). 

The integral f|T5|) is standard when A is positive semidefinite: 



y 

y 


r + 1' J 



(16) 



B 



(17) 



detA 



(19) 



Therefore, making the necessary algebraic simplifications and noting B^A ^B = — 4C, 
the purity of a coherent state is found to be 

This result shows that the entanglement for a two-mode coherent state in the center- 
of- mass/relative coordinates does not depend on the complex displacements a and /3, 
as expected. Since the displacement operators D{a,a) and D{f3,b) are separable in the 
atomic coordinates, then every coherent state must have the same entanglement as the 
ground state |0, 0). As a consequence, the purity of the ground state Pqo is also given 



by m 
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Figure 3. The ground state/coherent state entanglement Pqo as a function of /xi for 
four different values oi g = oj/V,: g = 1 (solid), g = 10 (dashed), g = 100 (dotted), and 
g = 1000 (dot-dashed). 



Figure 3 depicts the coherent state purity as functions of the kinematic parameter 
A^i (/^2 = 1 ~ /^i) Ei'iid dynamical parameter g = u/Q: 



Pai3{g,tJ'i) = PQo{g,fii 



9 



(21) 



When the trap and molecular frequencies are the same {g = 1), there is no atomic 
entanglement for any value of the mass ratios, as one might expect from the separability 
of the Hamiltonian. Note that either transformation g — )■ g^^ or — )■ fi2 leaves Pqo 
unchanged. As g departs from one, the purity decreases and the entanglement grows 
without limit. Maximum entanglement for a given g occurs when the masses are equal 
/^i = f^2 = 1/2. When the mass ratios of the two atoms are far off balance, the 
entanglement decreases. 

For the case of an unbound molecule, the purity of the time- dependent bound 
molecule $(t)on with center-of-mass momentum uncertainty at t = and relative 
vibrational state m = can be calculated the same way as the bound case, except the 
matrix A (fT6|) is now time dependent and takes the form 



A(t) 









F V? + 7^ 



v 



z 



z 

.V* 



z 

F Vi + 7' 




z 

z* 


■^2, ,2 I ^,2 



(22) 



where z = -7V2 + l/2e*'^FVi/^2 and 
steps as above, we find that 

Pouit) 



7F 



F + 7 V 
tan^^ r. Using this and following the same 

(23) 



Note that Pouit) has its maximum value at the moment of minimum uncertainty 
t = Mt/ (^F^) = and then decreases, meaning entanglement increases as time evolves. 
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Figure 4. The entanglement Pou{t) of $(t)oM at time i = as a function of /ii for 
four different values of c = 7/r: c = 1 (solid), c = 3 (dashed), c = 10 (dotted), and 
c = 30 (dot-dashed). 



An interesting relation to note is that when 7 = ry7^iA*2; there is no entanglement 
between the atoms in the free molecule at time t — 0. Figure 4 depicts the entanglement 
at t — using the parameterization c — F/j instead oi g — uj/ ft to highlight these 
features. 

To conclude this section, we compare these results for Gaussian states to results 
from the covariance matrix approach. Gaussian states are fully characterized by their 
first and second moments, but as wc have shown, the first moments have no influence on 
entanglement properties since they can be removed by a local unitary transformation. 
Defining the vector of operators 

/ V2{X, - (Xi)) \ 



V2 



(Pi - (Pi)) 



(A)) / 



(24) 



V 



Mi 



1 











the elements of the covariance matrix can be calculated from = {RkRi + RiRk)/'^- 
For either the ground state or the coherent state of the bound molecule, the covariance 
matrix evaluates to 

..2 



7^ + r V? 


A symplectic transformation to a new covariance matrix V = SVS^ with S G Sp(4, R) 
of the form S = S2 © S2 (S2 G Sp(2,]R)) is local with respect to the atomic observables 
and will not change the entanglement properties of the covariance matrix. In particular, 



r2 



V 







1 I 
+ ^2 





r2 



\ 

■72 + r ViAi2 


7^ + rv^ 



(25) 
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one can define a symplectic transformation that rescales the variables as 

/tTs \ ^ I TtT/s 



I/VtTs J ^ \ s/v^ ' ^'^^^ 

with = (7^ + r^/x^)/(7^ + r^/i2). This transformation brings the covariance matrix 
V into the standard form for a two-mode squeezed state [B] 



V 



/ cosh r 


sinhr 

V 





coshr 


— sinh r 



sinhr 


coshr 




\ 

— sinh r 


cosh r J 



(27) 



where the squeezing parameter r is directly related to the purity Pqo by cosh r = (Pqo)"^- 
For comparison with other results, note that for the covariance matrix V one finds 
that the logarithmic negativity Sj^ (a standard measure of entanglement for Gaussian 
states [5]) is exactly the squeezing parameter S^f = r. A similar result (but with a more 
complicated, time-dependent symplectic transformation S) holds for the unbound state 



4. Entanglement in Number States and Superpositions of Number States 

The purity for number states can be calculated using the connection between coherent 
states and number states 

1 (9" 



\n) 



da' 



-e 2 



la) 



(2^ 



a=0 



By substituting the wave functions for coherent state in atomic coordinates into the 
expression (128|) . we find 



■|/3|" 



From here we can proceed in two ways. We can substitute 
purity expression 



.(29) 

a,/3=0 

directly into the 



); that will be done below as part of the main text. An alternate 
approach uses an expansion onto another double harmonic oscillator basis, one that is 
separable in atomic coordinates and denoted |j, k}\ 

{Xi,X2\j,k} = {xi\j}{x2\k} 

1/2 



7i72 



Hj{'yiXi)Hk{'y2X2)e 



(30) 



2J+^7rj!A;!^ 

In this expression, 71 and 72 do not have a dynamical meaning based on the Hamiltonian 
like 7 and F, but instead are free parameters that should cancel out in the final expression 
for the purity. Applying f l28|) to the atomic basis vectors |j, A;}, we can find expressions 
for the coefficients (n, m |j. A;} that transform between the molecular number basis to 
the (artificial) atomic number basis. More details on this approach, which may be more 
useful for numerical simulations and for calculating the entropy of entanglement or other 
entanglement measures, are located in Appendix A. 
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into the purity expression 
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and 

(31) 



The real symmetric matrix A is the same as f|T6|) above, but now 

f 7(ai + «2) + r/ii(/3i + /32) \ 
7(a3 + «4) + ryUi(/33 + (3i) 
-7(ai + 04) + r/i2(/3i + /34) 
\ -7(03 + as) + r/i2(/33 + /32) / 

and 



B = 



(32) 



C = -l/2{al + al + al + al + + pl + /J^ + /?. 
Simplification leads to 



Pr 



00 



n 



z^Mz 



I). (33) 



(34) 



{a,,/3i}=0 



where 

z"^ = («!, 02, "3, "4, /32, /^S, /34) 

and M is an 8 X 8 matrix ( IB.3P described in Appendix B. 



(35) 



The expression f l34p can be evaluated analytically for all values of m and n, but it 
grows in complexity rapidly. Explicit calculations reveal that all Pmn have the form 

2m+2n 

E^{mn) 2ip4m+4n-2i (oc\ 



DO 



[(7^ + rv)(7' + rV2)]^ 

where C^"^"'^ are polynomials of /xi and /i2 with rational coefficients that can be 
determined from M. As examples, one finds 

7r 



01 



ry^ X (37^ + A^^'T^fll + III) 

4[(7^+/i?n)(72+^2r2)]5/2 W 

+ 2^'V\2^^\ + + 2/x^) + A^^V' + /i^) + ?>V'^^\^^1^ (37) 



and 



Pii 



7r 



16[(72 + /i?r2)(72 + ;,2r2)]9/2 

+ 127i2r^(8/it _ 3^2^2 ^ g^4^ ^ 2407'°rV?/i2(/i? + t^l) 
+27*^r*^(8/i? - 64/i?/i2 + 459^4^4 _ 64/i2/i6 + 8/i^) 
+2407'riVf/"2(/"? + lA) + 127'riVf/i2(8/ii - ^tAlA + 8/i^) 

+i6f riV?/i^(^? + tA) + sr^X^^) . (38) 

As for coherent states, one can rewrite any P^n so it depends only on the ratio of the 
frequencies g (or equivalently momentum uncertainty ratio c) and the mass fraction /ii 
(or ^2)- The lowest nine combinations of m, n are depicted in Figure 5. 
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Figure 5. Each column displays the graphs of Pmn for m = {0,1,2} (left, middle, 
right) and n — {0, 1, 2, 3} (solid, dashed, dotted, dot-dashed) for two values oi g: g = 1 
(top) and g — 5 (bottom). 



The following properties of the functions Pmn can be inferred either by analytical 
means or by graphical inspection of the calculated purity for the lowest combinations 
of m, n: 

• The purity functions are symmetric under exchange of trap and molecular quantum 
number, i.e. Pmn = Pnm- Similar to the ground state purity Pqq, the 
transformations /ii — )■ /i2 or (7 — )■ leave Pmn invariant, as can be expected 
from the symmetries of the Hamiltonian. 

• For all {m, n}, the energy eigenstates have no entanglement in the limiting cases of 
extreme mass difference, either m2 ^ mi (/ii ^ 0) or m2 ^ mi [fii ^1). In this 
unphysical limit, the length scale 7^^ = ^/h/Jiuj diverges and the energy eigenstates 
^m,n{xi,X2) become unnormalizable. 

• For all finite mass ratios, one finds Pmn < 1 and therefore the energy eigenstates are 
entangled. The only exception is the special case g = 1 and then only the ground 
state is separable. Inspecting the Hamiltonian in atomic coordinates ([3]) one can 
see that operator becomes separable whenever g = 1. However, this does not imply 
that the energy eigenfunctions $m,n(a;i, X2) become separable in that limit. As a 
side note, this does mean that for uncoupled oscillators there exists an alternate 
energy eigenstate basis \m, n) constructed of entangled states that coincides with 
the atomic oscillator basis |j, k} only on the ground state. 

• When (7 = 1, Pmn is a polynomial of /ii of order 2(m + n). For example, Pqo = 1) 
P^o = Poi = 1 - 2/ii + 2/i2, and = 1 - 8/ii + 32/i2 - 48/i? + 24/i^. 

• Based on graphical analysis, for a fixed value of g it appears Pm+i,n+i < Pm.,n for all 
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/ii, but no other inequalities appear universal. For example, there are some regions 
of {(7, /ii}-parameter space where P12 > Pu and where P13 > P12. As m and n 
increase the functions become more oscillatory and more tightly spaced, making it 
unlikely to hypothesize any other bounds based on graphical methods alone. 

• For a fixed value of fii, the entanglement generally increases as g increases. As 
(? — i- 00, the purity takes its global minimum at /ii = 1/2, although for general 
g, Pmn need not have a minimum at fii = 1/2, and it may in fact have a local 
maximum. For fixed /ii, the entanglement increase is not monotonic in g, but 
shows local maxima and minima of decreasing prominence as g increases. 

As before, one can use the same method, only replacing A (fT6!) with A(t) (!22|) . 
to calculate $(t)m« for unbound diatomic molecules with center-of-mass momentum 
uncertainty hT/\/2 and vibrational state m. 

Finally, in principle one can calculate the entanglement of general states 

00 

1"^) = Cm,n\m,n) (39) 



m,n=0 



as 



where 



{mi,ni}ig{i_2,3,4}=0 



P{{m„n,}) = ^^ I n e^"^^ (41) 



In contrast to the purities of number states, Pmn, the functions P({mj,?2j}) are not 
necessarily positive or symmetric around fii = 1/2. Nor do they limit to unity when 
/ii — 7- or 1, although numeric analysis suggests they are bounded functions. Many 
symmetry relations between permutations of indices can be derived. For example, 
because of the integration range, only even kernels contribute, meaning P({mj, rii}) = 

unless J2t=ii^i + ^i) is even. 

Figure 6 gives a flavor for the entanglement properties of molecular state super- 
positions. The purity of the state 

1$) = cose|0,l) + sine|l,0) (42) 

is depicted as a function of fii for several values of g and 6. 



5. Conclusions: Correlation, Covariance, and the Classical-Quantum 
Correspondence 

When the trap frequency and molecular frequency are the same {g = 1), there is no 
entanglement between the atoms when the center-of-mass and relative oscillators are 
in coherent states, and in particular the ground state. For all other coherent states 
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Figure 6. This figure depicts tlie purity of tlie reduced density matrix of for two 
values of g and three values of 6*: 9 = ?i7r/2 (solid), 9 = mr/2 + tt/6 (dashed), and 
9 = mr/2 + 7r/3 (dashed) for n any integer. Note that for g — I and 9 ^ mT/2, there 
are certain mass ratios that have no interatomic entanglement. 



and number states with finite mass ratios, there is entanglement between the atoms, 
which generally increases with g and with the index of the number states. Certain 
linear combinations of number states, including those entangled with respect to the 
molecular observables, can be disentangled with respect to atomic observables, but 
again apparently only when g = 1. These mathematical results, based on the definition 
of separability with respect to a given tensor product structure, can be mathematically 
proven. 

In interacting systems, correlations are expected, and one perspective is that this 
kind of entanglement is just an artifact of studying the molecule in the 'wrong' basis, 
i.e. the atomic basis. Since the Hamiltonian is separable with respect to the center-of- 
mass/relative basis, in some sense the dynamics 'chooses' the molecular observables over 
the particle observables. However, if the atoms have internal structure, one can imagine 
at least in principle, that the atomic observables could be experimentally accessed (or 
"chosen") independently of the molecular observables and attempts could be made to 
measure correlations between the atoms. For example, one could attempt to measure 
the covariance between Xi and X2: 

(T.,x, = V13/2 = - (Xi)(X2). (43) 

For molecular coherent states, one finds 

1 /ii/ii 1 _i 



272 ^''^ 3 )■ '^^^^ 

This covariance quantifies the correlations in uncertainty that the atomic positions 
inherit from the intrinsic minimum uncertainty of the molecular oscillators. Low 
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frequency traps with high frequency molecules imply the largest covariance (and 
therefore easiest to measure). When g = 1, the covariance disappears for coherent 
states, as does the entanglement, which can be seen from ( 125|) . The same holds for 
the other elements of the covariance matrix. However, although there is entanglement, 
one can show that these correlations are not intrinsically quantum. We could imagine 
a classical analogue: two masses on a spring. The masses could be at rest, but with 
statistical uncertainty in their positions, so that the probability distribution as a function 
of center of mass coordinates {x,p} and relative coordinates {r,q} is 

p(x,p,r,q) = _l-e-r^(---o)^e-^^^"^°^'e-^'(^-^°)'e-^(^-''°^'. (45) 

This supposition gives the same covariance matrix V as the quantum coherent state. 
The only difference between the quantum and classical covariance is that in the quantum 
case 7 and F are of dynamical origin, whereas they are of purely statistical origin in the 
classical case. Another way of saying this is that the transformation from molecular 
observables to atomic observables maps a positive, Gaussian Wigner function into 
another positive Gaussian Wigner function, and one can show that some classical model 
can reproduce correlations in any positive Wigner function. 

In contrast, in the number state |m,n) the position covariance (Tx-i^x2 is 

^.... = ^(2r^ + l)-|^(2m + l) 

= ^((2n + l)-r'(2m + l)). (46) 

This has very different properties from the coherent state covariance, and from the 
entanglement of such states. When n = m and for g = 1 the covariance in number 
states disappears, although the entanglement is non-vanishing for all m and n except 
m = n = 0. The covariance is greatest in magnitude for (7 ^ 1 and large n or 
g 1 and large m, a relationship between g and number that does not exist for 
the entanglement. For number states there is no correspondence to a classical model 
with statistical uncertainty. Additionally, we note that the Wigner function for a number 
state has negative regions, often considered a signal of 'quantumness', whether expressed 
in atomic or molecular variables, but entanglement only in atomic variables. 

As a final comment, we note that the expressions for purity do not depend on h or 
the ratio of the quantum scale to the classical scale in any way. Although measurements 
of the postion and momentum uncertainties will involve a scale set by h, in the purity 
expression only the ratios appear and so the overall scale cancels out. To see whether 
the entanglement correlations are truly "quantum" one could also imagine constructing 
dichotomous observables on the (xi, X2)-space such that Bell-type inequalities can be 
formulated. Such observables can be constructed in several ways, for example based 
on the displaced parity operator [23] or using pseudospin operators [25], and this 
investigation will be pursued in future work. 
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Appendix A. An Alternate Approach 

If the goal is to calculate the entanglement of a state with respect to the atomic tensor 
product structure "H = "Hi 7/2, then any basis that is separable with respect to 
this structure can be used for taking the partial trace. In the main body of this 
paper, the purity of the reduced density matrices was calculated using the continuous- 
variable atomic coordinate basis |a;i,X2}. In this appendix, we instead use a double 
harmonic oscillator basis |j, /c}. These states are realized by separable wave functions 
( l30l) characterized by positive real parameters 71 and 72 that can be freely chosen 
for convenience. There can be advantages of using such a discrete basis to generate 
approximate expressions for the purity even when an exact analytic expression can also 
be derived. 

Using the |j, k} basis, the reduced density matrix for atom 1 can be written 

00 

Pi= (A.i) 

and the purity of the reduced density matrix is 

00 

P($)= {jM^){^\fM{f^'\^){^\j,k'}. (A.2) 

j,j',k,k'=a 

We will focus on calculating the entanglement for number states |$) = |m, n), 
and the matrix elements transforming between the molecular oscillators and the atomic 
oscillators {j, k\m,n) are the central objects of concern. Similar to the procedure in the 
main text, an expression for this matrix element will be derived by taking the derivatives 
of the matrix element between coherent states {Ti,r2| a, (3): 

{j, k\ m, n) = ] (A.3) 



Qj+k+m+n 



ri,T2,a,/3=0 



dTidT^da"'d(5' 

The molecular coherent state |a,/3) is defined above ([7]) and the atomic coherent states 
are defined as 

00 j 

|n,r2} = e-l-l^/-l-l^/^5:^b,A;}. (A.4) 

j,k=0 ^ 

The associated wave function in particle coordinates is 
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X e-^(2:i-na:)^-^{a;2-T2:r)^g|{ripXl+T2pX2)^ (A. 5) 

with analogs definitions for ti^, Tip, etc., to (fl^ . One can tlien evaluate {ri, a, (3) by 
performing the integral 

{ti, 7-2I a, /3) = / rfa;i(ix2 {ri,r2|xi,a;2)(xi,a;2| (A.6) 



This is yet another Gaussian integral, but this time only in two variables. Completing 
the integration and substituting into (1A.3I) . the coefficient can be written 

1/2 

/ i-i-' u -I ■ V n' u I \ 

{i,k\m,n) 



/ 47i727r 
\^j\k\m\n\Z 



aj+k+m+n 

X -Jl g-l/2(r2+r|+a2+/32)gF/Z 



(A.7) 



'''l,T2,a,l3=0 



where 



F = (72 + /i^r^ + 72)(7a + /iir/3 + 7iri)2 

+ 2(7^ - /ii/i2r^)(7" + Atir/3 + 7in)(-7" + 1^2^ P + 12T2) 

+ (7^ + filT^ + 7?)(-7« + /i2r/3 + 721-2)^ 

Z = 7'r2 + /,2^2r2 ^ ^2^2p2 ^ ^2^2 ^ ^2 ^^2 ^ ^2)_ (^^g) 

The expression (1A.7P depends in a complicated fashion on the non-physically 
meaningful parameters 71 and 72, but surprisingly, when this coefficient is substituted 
into the summation in flA.2p . this dependence must cancel. To examine how this sum 



converges to the exact result flMl) . Figure 7 depicts the convergence of Pqi ioi four 
cases of {g,fii) and for four several different values of (71,72). We do not answer the 
potentially interesting question of how to choose 71 and 72 for optimal convergence, 
although some features make intuitive sense. For example, for smaller g, smaller choices 
for (71,72) converge faster. Also, for the cases where /ii = 1/10, choices with 71 < 72 
converge faster than those with 71 > 72. 

The advantage of using this method is that the reduced density matrix pi can 
be approximated to arbitrary accuracy by a finite-dimensional matrix. This matrix, for 
example, could be diagonalized and used to calculate the entropy of entanglement. Also, 
the coefficients {j, k\m,n) are necessary if one is to use the pseudo-spin operators [25] 
to construct Bell-type inequalities for this system. More generally, even though one 
can find exact analytic continuous-variable expressions for the purity of the reduced 
matrix elements, numerical schemes require discretization, which in certain cases has 
been shown to mask the presence of entanglement in continuous variable systems [T8] . 



Appendix B. The matrix M 

The entanglement of an number state can be calculated exactly using the expression 



Poo ^ ^,TM, 



(B.l) 

{a^,|3i}=0 
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Figure Al. This figure depicts the convergence of the approximation scheme for 
Poi for four different values of {g,Hi) and four different choices of the parameters 
(71:72): (71,72) = (l/\/2, 1/V2) (circle), (1,1) (square), (1/^2,1) (up-triangle), and 
(l,l/-\/2) (down-triangle). For each combination of (g,/ui) and (71,72), the accuracy 
of six successive approximations to Pqi are plotted. Each approximation corresponds 
to taking more and more terms in the sum in (|A.2p from jmax = fcmax = to 5. If not 
all six shapes are apparent, successive approximation give results indistinguishable on 
this scale. For the first case, {g,fJ.i) = (1, 1/2), the first choice (71,72) — (l/\/2, l/"\/2) 
gives the exact result at the second approximation jmax = ^max — 1 • 



where 

= (ai, a2, as, a^, /32, /^s, f^i)- 
The matrix M is 8 x 8 and can be written as 
/ 



M 



where 
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s = 7r(7^ - r^/ii/i2)(Aii - Ai2)/-D 



(B.2) 



(B.3) 



71(7^ + T^fiifi2)/D 
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L> = 4(72 + rVi)(7' + rV2) 

It is useful to note that det M — 1/256. 



